/*=============================================================================
*
*   Copyright (C) 2020 All rights reserved.
*
*   Filename:		levelset.hpp
*
*   Author: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Date: 2020-03-23 04:16
*
*   Last Editors: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Last modified:	2020-08-04 01:25
*
*   Description:
*
=============================================================================*/
#ifndef LEVELSET_H
#define LEVELSET_H

#include "../arithmetic/basic_arithmetic.hpp"
#include "../types.hpp"
#include <boost/geometry.hpp>

namespace geoxfem {
namespace xfem {

template <typename point_type, typename CrackTip>
struct get_level_set { /*{{{*/
  static std::pair<_real, _real> apply(const point_type& pt,
                                       const CrackTip&   tip) {
    _real phi, psi;

    /* if point lies in front of tip */
    VectorXR tangent_dir = tip.getTangentVector();
    VectorXR normal_dir  = tip.getNormalVector();
    VectorXR point_vec   = pt.getCoord() - tip.getTip().getCoord();
    psi                  = tangent_dir.dot(point_vec);
    if(psi > 0) {
      phi = tip.getNormalVector().dot(point_vec) * tip.getSign();
      return std::make_tuple(psi, phi);
    }

    /* if point lies behind the tip */
    psi = 0;
    std::vector<std::pair<_real, _real>> candidates;
    for(auto pc0 = tip.begin(), pc1 = std::next(pc0); pc1 != tip.end();
        ++pc0, ++pc1) {
      tangent_dir = (pc1->getCoord() - pc0->getCoord()).normalized();
      normal_dir << -tangent_dir[1], tangent_dir[0];
      point_vec        = pt.getCoord() - pc0->getCoord();
      int   side       = arithmetic::sgn(point_vec.dot(normal_dir));
      _real projection = std::abs(point_vec.dot(tangent_dir));
      _real length     = boost::geometry::distance(*pc0, *pc1);
      phi              = side
            * boost::geometry::distance(pt, std::vector<point_type>{*pc0, *pc1})
            * tip.getSign();
      candidates.emplace_back( std::make_pair(psi - std::min({projection, length}), phi)); psi -= length;
       }
    auto it = std::min_element(
      candidates.begin(), candidates.end(),
      [](const std::pair<_real, _real>& v1, const std::pair<_real, _real>& v2) {
        return std::abs(v1.second) < std::abs(v2.second);
      });
    return *it;
  }
}; /*}}}*/

template <typename point_type, typename CrackTip>
struct get_level_set<point_type, CrackTip*> { /*{{{*/
  static std::pair<_real, _real> apply(const point_type& pt,
                                       const CrackTip*   ptr_tip) {
    return get_level_set<point_type, CrackTip>::apply(pt, *ptr_tip);
  }
}; /*}}}*/

template <typename PointLike, typename CrackTip>
std::pair<_real, _real> GetLevelSet(const PointLike& pt,
                                    const CrackTip&  crack) {
  return get_level_set<PointLike, CrackTip>::apply(pt, crack);
}

}  // namespace xfem
}  // namespace geoxfem

#endif /* LEVELSET_H */
